Priors
\(\alpha\)
tibble(x = seq(0, 1.3, length = 10^(5)),
y = with(prior_params, gamma_density(x,
mean = alpha_mean,
sd = alpha_sd,
bounds = alpha_bounds))) %>%
ggplot(aes(x=x, y = y)) +
geom_line(alpha = .8) +
geom_ribbon(aes(x=x,ymin=0,ymax=y),
fill="black",
alpha=.6) +
theme_c(legend.text=element_text(size = 10)) +
viridis::scale_fill_viridis(discrete=TRUE,
option = "rocket", begin=.3,end=.8) +
labs(x = "Value",
y = "Probability Density",
title = latex2exp::TeX("Definition of Prior for $\\alpha$"),
fill ='',
subtitle = paste0("Mean: ", prior_params$alpha_mean,
", SD: ", prior_params$alpha_sd))
\(\beta\)
tibble(x = seq(0, 1, length = 10^(5)),
y = with(prior_params, beta_density(x,
mean = beta_mean,
sd = beta_sd,
bounds = beta_bounds))) %>%
ggplot(aes(x=x, y = y)) +
geom_line(alpha = .8) +
geom_ribbon(aes(x=x,ymin=0,ymax=y),
fill="black",
alpha=.6) +
theme_c(legend.text=element_text(size = 10)) +
viridis::scale_fill_viridis(discrete=TRUE,
option = "rocket", begin=.3,end=.8) +
labs(x = "Value",
y = "Probability Density",
title = latex2exp::TeX("Definition of Prior for $\\beta$"),
fill ='',
subtitle = paste0("Mean: ", prior_params$beta_mean,
", SD: ", prior_params$beta_sd))
\(P(S_1|\text{untested})\)
tibble(x = seq(0, 1, length = 10^(5)),
y = with(prior_params, beta_density(x,
mean = s_untested_mean,
sd = s_untested_sd,
bounds = s_untested_bounds))) %>%
ggplot(aes(x=x, y = y)) +
geom_line(alpha = .8) +
geom_ribbon(aes(x=x,ymin=0,ymax=y),
fill="black",
alpha=.6) +
theme_c(legend.text=element_text(size = 10)) +
viridis::scale_fill_viridis(discrete=TRUE,
option = "rocket", begin=.3,end=.8) +
labs(x = "Value",
y = "Probability Density",
title = latex2exp::TeX("Definition of Prior for $P(S_1|untested)$"),
fill ='',
subtitle = paste0("Mean: ", prior_params$s_untested_mean,
", SD: ", prior_params$s_untested_sd))
tibble(x = seq(0, 1, length = 10^(5)),
y = with(prior_params, beta_density(x,
mean = s_untested_mean,
sd = s_untested_sd,
bounds = s_untested_bounds))) %>%
ggplot(aes(x=x, y = y)) +
geom_line(alpha = .8) +
geom_ribbon(aes(x=x,ymin=0,ymax=y),
fill="black",
alpha=.6) +
theme_c(legend.text=element_text(size = 10)) +
viridis::scale_fill_viridis(discrete=TRUE,
option = "rocket", begin=.3,end=.8) +
labs(x = "Value",
y = "Probability Density",
title = latex2exp::TeX("Definition of Prior for $P(S_1|untested)$"),
fill ='',
subtitle = paste0("Mean: ", prior_params$s_untested_mean,
", SD: ", prior_params$s_untested_sd))
\(P(S_0| \text{test}_+, \text{untested})\)
tibble(x = seq(0, 1, length = 10^(5)),
y = with(prior_params, beta_density(x,
mean = p_s0_pos_mean,
sd = p_s0_pos_sd,
bounds = p_s0_pos_bounds))) %>%
ggplot(aes(x=x, y = y)) +
geom_line(alpha = .8) +
geom_ribbon(aes(x=x,ymin=0,ymax=y),
fill="black",
alpha=.6) +
theme_c(legend.text=element_text(size = 10)) +
viridis::scale_fill_viridis(discrete=TRUE,
option = "rocket", begin=.3,end=.8) +
labs(x = "Value",
y = "Probability Density",
title = latex2exp::TeX("Definition of Prior for $P(S_0|test_+,untested)$"),
fill ='',
subtitle = paste0("Mean: ", prior_params$p_s0_pos_mean,
", SD: ", prior_params$p_s0_pos_sd))
tibble(x = seq(0, 1, length = 10^(5)),
y = with(prior_params, beta_density(x,
mean = .55,
sd = .12,
bounds = NA))) %>%
ggplot(aes(x=x, y = y)) +
geom_line(alpha = .8) +
geom_ribbon(aes(x=x,ymin=0,ymax=y),
fill="black",
alpha=.6) +
theme_c(legend.text=element_text(size = 10)) +
viridis::scale_fill_viridis(discrete=TRUE,
option = "rocket", begin=.3,end=.8) +
labs(x = "Value",
y = "Probability Density",
title = latex2exp::TeX("Definition of Prior for $P(S_0|test_+,untested)$"),
fill ='',
subtitle = paste0("Mean: ", .55,
", SD: ", .12)) +
scale_x_continuous(n.breaks=11)
qbeta(.025, shape1= get_beta_params( mu = .55,
sd = .12)$a,
shape2= get_beta_params( mu = .55,
sd = .12)$b
)## [1] 0.3124573qbeta(.975, shape1= get_beta_params( mu = .55,
sd = .12)$a,
shape2= get_beta_params( mu = .55,
sd = .12)$b
)## [1] 0.7758402COVID-19 Trends and Impact Survey
\(\beta\)
ctis_smoothed %>%
filter(keep) %>%
mutate(state=toupper(state)) %>%
select(date,
state,
imputed_beta,
beta_estimate_smoothed,
beta_estimate_spline_smoothed) %>%
pivot_longer(contains("beta")) %>%
mutate(name = case_when(
name == "beta_estimate_smoothed" ~ "LOESS smoothed",
name == "beta_estimate_spline_smoothed" ~ "Spline smoothed",
name == "imputed_beta" ~ "Survey Value"
)) %>%
ggplot(aes(x=date, y=value, color = name,
alpha = name, linewidth=name)) +
geom_line() +
facet_wrap(~state, ncol=4) +
scale_alpha_manual(values=c("LOESS smoothed" = .9,
"Spline smoothed" = .9,
"Survey Value"=.3),
name='') +
scale_linewidth_manual(values = c("LOESS smoothed" = 1.05,
"Spline smoothed" = 1.05,
"Survey Value"=.5)) +
scale_color_manual(values=c("#3381FF", "#B58746", "#26900F"),
name ='') +
guides(alpha="none",
linewidth="none",
color = guide_legend(override.aes = list(linewidth = 3,
alpha = c("LOESS smoothed" = .9,
"Spline smoothed" = .9,
"Survey Value"=.3)),
nrow=3)) +
theme_c(legend.position="top") +
ylim(0,1) +
scale_x_date(date_breaks="3 months",
date_labels = "%b %Y") +
labs(y = TeX("Survey Estimate of $\\beta$"),
x= "",
title = TeX("Comparing Approaches for Smoothing Survey Estimates of $\\beta$"))
The ratio of the screening test positivity over the overall test positivity from the COVID-19 Trends and Impact Survey is taken to be the estimate of \(\beta\). We compare two approaches to smoothing: cubic spline smoothing with 2 knots (July 15th, 2021 and December 1st, 2021) to LOESS smoothing with a span of 0.33.
num_states <- ctis_smoothed %>%
filter(keep) %>%
pull(state) %>%
unique() %>%
length()
cat("There are", num_states, "states with sufficient data available to estimate $\\beta$ across the time period considered.")## There are 27 states with sufficient data available to estimate $\beta$ across the time period considered.ctis_smoothed %>%
filter(keep) %>%
mutate(state=toupper(state)) %>%
select(date,
state,
imputed_beta,
# beta_estimate_smoothed,
beta_estimate_spline_smoothed) %>%
pivot_longer(contains("beta")) %>%
mutate(name = case_when(
name == "beta_estimate_spline_smoothed" ~ "Spline smoothed",
name == "imputed_beta" ~ "Survey Value"
)) %>%
ggplot(aes(x=date, y=value, color = name,
alpha = name, linewidth=name)) +
geom_line() +
facet_wrap(~state, ncol=4) +
scale_alpha_manual(values=c(
"Spline smoothed" = .9,
"Survey Value"=.3),
name='') +
scale_linewidth_manual(values = c(
"Spline smoothed" = 1.05,
"Survey Value"=.5)) +
scale_color_manual(values=c("#3381FF", "#26900F"),
name ='') +
guides(alpha="none",
linewidth="none",
color = guide_legend(override.aes = list(linewidth = 3,
alpha = c(
"Spline smoothed" = .9,
"Survey Value"=.3)),
nrow=3)) +
theme_c(legend.position="top") +
ylim(0,1) +
scale_x_date(date_breaks="3 months",
date_labels = "%b %Y") +
labs(y = TeX("Survey Estimate of $\\beta$"),
x= "",
title = TeX("Smoothing Survey Estimates of $\\beta$"))
The ratio of the screening test positivity over the overall test positivity from the COVID-19 Trends and Impact Survey is taken to be the estimate of \(\beta\). Because of the noise present inthe data, we use the spline smoothed values to inform the priors. More specifically, we use cubic spline smoothing with 2 knots (July 15th, 2021 and December 1st, 2021).
\(\Pr(S_1 | \text{untested})\)
ctis_smoothed %>%
mutate(state=toupper(state)) %>%
select(date,
state,
contains("s_untested")) %>%
pivot_longer(contains("s_untested")) %>%
mutate(name = case_when(
name == "s_untested_smoothed" ~ "LOESS smoothed",
name == "imputed_s_untested" ~ "Survey Value"
)) %>%
ggplot(aes(x=date, y=value, color = name,
alpha = name, linewidth=name)) +
geom_line() +
facet_wrap(~state, ncol=4) +
scale_alpha_manual(values=c("LOESS smoothed" = .9,
"Survey Value"=.6),
name='') +
scale_linewidth_manual(values = c("LOESS smoothed" = 1.02,
"Survey Value"=.9)) +
scale_color_manual(values=c("Survey Value" = "#26900F",
"LOESS smoothed"="#3381FF"),
name ='') +
guides(alpha="none",
linewidth="none",
color = guide_legend(override.aes = list(
linewidth = 3,
alpha = c("LOESS smoothed" = .9,
"Survey Value"=.3)),
nrow=3)) +
theme_c(legend.position="top") +
scale_x_date(date_breaks="3 months",
date_labels = "%b %Y") +
labs(y = TeX("Survey Estimate of $Pr(S_1|untested)$"),
x= "",
title = TeX("Smoothing Survey Estimates of $Pr(S_1|untested)$"))
The percentage of the population experiencing COVID-19-like illness is taken to be the estimate of $(S_1|). The LOESS smoothed estimate with a span of 0.2 is shown in red.
Distributions of \(\beta\) by State
# raw ctis data
ctis_raw <- readRDS(here('data/data_raw/ctis_all_states.RDS'))
states_keep <- ctis_smoothed %>%
filter(keep) %>% pull(state) %>% unique()
ctis <- ctis_raw %>%
select(signal, state=geo_value, date, value) %>%
pivot_wider(names_from=signal, values_from=value) %>%
filter(state %in% states_keep) %>%
mutate(beta=smoothed_wscreening_tested_positive_14d/
smoothed_wtested_positive_14d,
state=toupper(state))ctis %>%
group_by(state) %>%
summarize(median = median(beta, na.rm=TRUE),
mean = mean(beta, na.rm=TRUE),
q1 = quantile(beta, .025, na.rm=TRUE),
q2 = quantile(beta, .975, na.rm=TRUE)) %>%
ggplot(aes(x= fct_reorder(state, mean, .desc=TRUE),
y = mean)) +
geom_errorbar(aes(ymin=q1,ymax=q2),width=.2) +
geom_point(color="darkred", size=2) +
scale_y_continuous(n.breaks=10) +
labs(y= TeX("Empirical Estimates of $\\beta$"),
x="State",
title = TeX("Empirical Estimates of $\\beta$ from the COVID-19 Trends and Impact Survey"),
subtitle = "2.7% Percentile, 97.5% Percentile, and Median By State" ) +
theme_c(axis.title=element_text(size=12))
Since \(\beta\) represents the ratio of \(\Pr(\text{test}_+|\text{untested}, S_0)\) to \(\Pr(\text{test}_+|\text{tested})\), we estimate \(\beta\) from the COVID-19 Trends and Impact Survey as the ratio of the screening test positivity over the overall test positivity. Here, we consider the 2.5% and 97.5% percentiles and mean for each state. Most have a mean near 0.2.
#
# Logistic transformation of the Beta CDF.
#
f.beta <- function(alpha, beta, x, lower=0, upper=1) {
p <- pbeta(x, alpha, beta)
log(p/(1-p))
}
#
# Sums of squares.
#
delta <- function(fit, actual) sum((fit-actual)^2)
#
# The objective function handles the transformed parameters `theta` and
# uses `f.beta` and `delta` to fit the values and measure their discrepancies.
#
objective <- function(theta, x, prob, ...) {
ab <- exp(theta) # Parameters are the *logs* of alpha and beta
fit <- f.beta(ab[1], ab[2], x, ...)
return (delta(fit, prob))
}
#
# Solve two problems.
#
par(mfrow=c(1,2))
alpha <- 15; beta <- 22 # The true parameters
for (x in list(c(1e-3, 2e-3), c(1/3, 2/3))) {
x.p <- f.beta(alpha, beta, x) # The correct values of the p_i
start <- log(c(1e1, 1e1)) # A good guess is useful here
sol <- nlm(objective, start, x=x, prob=x.p, lower=0, upper=1,
typsize=c(1,1), fscale=1e-12, gradtol=1e-12)
parms <- exp(sol$estimate) # Estimates of alpha and beta
#
# Display the actual and estimated values.
#
print(rbind(Actual=c(alpha=alpha, beta=beta), Fit=parms))
#
# Plot the true and estimated CDFs.
#
curve(pbeta(x, alpha, beta), 0, 1, n=1001, lwd=2)
curve(pbeta(x, parms[1], parms[2]), n=1001, add=TRUE, col="Red")
points(x, pbeta(x, alpha, beta))
}## alpha beta
## Actual 15.00000 22.00000
## Fit 14.99983 21.99811## alpha beta
## Actual 15.00000 22.00000
## Fit 14.99926 21.99926
ctis %>%
group_by(state) %>%
summarize(median = median(beta, na.rm=TRUE),
mean = mean(beta, na.rm=TRUE),
q1 = quantile(beta, .025, na.rm=TRUE),
q2 = quantile(beta, .975, na.rm=TRUE)) %>%
ggplot(aes(y=fct_reorder(state, mean, .desc=TRUE))) +
ggridges::geom_density_ridges(
aes(x=beta,
y=fct_reorder(state, beta, .fun =median)),
data=ctis,
fill="#596281",
quantile_lines=TRUE,
quantiles = c(.025,.975),
quantile_fun=function(beta,...) c(
# quantile(beta,.025, na.rm=TRUE),
median(beta,na.rm=TRUE)
# quantile(beta,.975, na.rm=TRUE)
)) +
# geom_errorbar(aes( xmin=q1,xmax=q2),width=.2) +
scale_x_continuous(n.breaks=10, limits=c(0,1)) +
labs(x= TeX("Empirical Estimates of $\\beta$"),
y="State",
title = TeX("Empirical Estimates of $\\beta$ from the COVID-19 Trends and Impact Survey"),
subtitle = "Density Estimate of Distribution and Median by State" ) +
theme_c(axis.title=element_text(size=12),
axis.text.x=element_text(size=11)) 
ctis %>%
ggplot(aes(x = beta, fill=state)) +
geom_density(alpha=.6) +
scale_fill_viridis(option="mako", discrete=TRUE) +
theme_c()
## [1] 0.2160426ctis %>%
group_by(state) %>%
summarize(mean=mean(smoothed_wcli,na.rm=TRUE)) %>%
pull(mean) %>%
mean()## [1] 0.01602321Compare Priors
source(here("R/05-state-analysis.R"))
beta_shape = get_shape(ctis_smoothed,
option="beta",
quiet=FALSE)
tibble(x = seq(0, 1, length = 10^(4)),
original = with(prior_params, beta_density(x,
mean = beta_mean,
sd = beta_sd,
bounds = beta_bounds)),
ctis= dbeta(x, shape1=beta_shape[1],
shape2=beta_shape[2])) %>%
pivot_longer(c(original,ctis))%>%
ggplot(aes(x=x, y = value, fill=name)) +
geom_line(alpha = .8) +
geom_ribbon(aes(x=x,ymin=0, ymax=value,fill=name),
alpha=.6) +
theme_c(legend.text=element_text(size = 10)) +
viridis::scale_fill_viridis(discrete=TRUE,
option = "rocket", begin=.3,end=.8) +
labs(x = "Value",
y = "Probability Density",
title = latex2exp::TeX("Definition of Prior for $\\beta$"),
fill ='')
s_untested_shape = get_shape(ctis_smoothed,
option="s_untested",
quiet=FALSE)
tibble(x = seq(0, 1, length = 10^(4)),
original = with(prior_params, beta_density(x,
mean = s_untested_mean,
sd = s_untested_sd,
bounds = s_untested_bounds)),
ctis= dbeta(x, shape1=s_untested_shape[1],
shape2=s_untested_shape[2])) %>%
pivot_longer(c(original,ctis))%>%
ggplot(aes(x=x, y = value, fill=name)) +
geom_line(alpha = .8) +
geom_ribbon(aes(x=x,ymin=0, ymax=value,fill=name),
alpha=.6) +
theme_c(legend.text=element_text(size = 10)) +
viridis::scale_fill_viridis(discrete=TRUE,
option = "rocket",
begin=.3,
end=.8) +
labs(x = "Value",
y = "Probability Density",
title = latex2exp::TeX("Definition of Prior for $Pr(S_1|untested)$"),
fill ='')
ctis_smoothed <- tar_read(ctis_smoothed,
store=here("_targets"))
ctis_smoothed %>% left_join(dates) %>%
group_by(biweek,state) %>%
summarize(beta_est=mean(beta_est,na.rm=TRUE)) %>%
ungroup() %>%
summarize(q1 = quantile(beta_est,.025, na.rm=TRUE),
median = median(beta_est, na.rm=TRUE),
q2 = quantile(beta_est,.975, na.rm=TRUE),
sd=sd(beta_est,na.rm=TRUE))
ctis_smoothed %>%
summarize(q1 = quantile(
beta_estimate_spline_smoothed,.025, na.rm=TRUE),
median = median(beta_estimate_spline_smoothed, na.rm=TRUE),
q2 = quantile(beta_estimate_spline_smoothed,.975, na.rm=TRUE),
sd=sd(beta_estimate_spline_smoothed,na.rm=TRUE))Bayesian Melding
library(patchwork)
library(truncdist)
source(here('R/04-melding.R'))
params <- prior_params
params$beta_shape1 <- beta_shape[1]
params$beta_shape2 <- beta_shape[2]
params$s_untested_shape1 <- s_untested_shape[1]
params$s_untested_shape2 <- s_untested_shape[2]
params$direct_params <- TRUE
#
# prior_params <- list(
# alpha_mean = .95,
# alpha_sd = 0.08,
# alpha_bounds = NA,
# # alpha_bounds = c(.8,1),
# beta_mean = .15,
# beta_sd =.09,
# beta_bounds = NA,
# # beta_bounds = c(0.002, 0.4),
# s_untested_mean = .03,
# s_untested_sd = .0225,
# # s_untested_bounds = c(0.0018, Inf),
# s_untested_bounds = NA,
# p_s0_pos_mean = .4,
# p_s0_pos_sd = .1225,
# # p_s0_pos_bounds = NA,
# p_s0_pos_bounds = c(.25, .7),
# pre_nsamp = 1e4,
# post_nsamp = 1e4)
melded <- do.call(get_melded, prior_params)
theta_wide <- with(prior_params,
theta <- tibble(x= seq(0,1.4, length = 200),
alpha = gamma_density(x,
mean = alpha_mean,
sd = alpha_sd,
bounds = alpha_bounds),
# beta= beta_density(x,
# mean = beta_mean,
# sd = beta_sd,
# bounds = beta_bounds),
# P_S_untested = beta_density(x,
# mean = s_untested_mean,
# sd = s_untested_sd,
# bounds = s_untested_bounds),
P_S0_pos = beta_density(x,
mean = p_s0_pos_mean,
sd = p_s0_pos_sd,
bounds = p_s0_pos_bounds)
)) %>%
mutate(beta = dbeta(x, shape1= beta_shape[1], shape2=beta_shape[2]),
P_S_untested = dbeta(x, shape1=s_untested_shape[1],
shape2=s_untested_shape[2] ))
theta_samp<- with(prior_params,
theta <- tibble(alpha = sample_gamma_density(pre_nsamp,
mean = alpha_mean,
sd = alpha_sd,
bounds = alpha_bounds),
beta= rbeta(pre_nsamp,
shape1= beta_shape[1],
shape2=beta_shape[2]),
P_S_untested = rbeta(pre_nsamp,
shape1=s_untested_shape[1],
shape2=s_untested_shape[2])) %>%
mutate(phi_induced = est_P_A_testpos(P_S_untested = P_S_untested,
alpha = alpha,
beta=beta))) %>%
select(phi_induced)
theta <- theta_wide %>%
pivot_longer(-c(x)) %>%
mutate(name = case_when(
name == "alpha" ~"$\\alpha$",
name == "beta" ~"$\\beta$",
name == "phi_induced" ~ "$M(\\theta) = Pr(S_0|test+,untested)$",
name == "P_S_untested" ~ "$Pr(S_1|untested)$",
name == "P_S0_pos" ~ "$Pr(S_0|test+,untested)$")
) %>%
mutate(name = factor(name,
levels = c(
"$\\alpha$",
"$\\beta$",
"$Pr(S_1|untested)$",
"$Pr(S_0|test+,untested)$",
"$M(\\theta) = Pr(S_0|test+,untested)$"))) melded_long <- reformat_melded(melded$post_melding,
melded$pre_melding,
pre_nsamp = prior_params$pre_nsamp,
p_s0_pos_mean = prior_params$p_s0_pos_mean,
p_s0_pos_sd = prior_params$p_s0_pos_sd,
p_s0_pos_bounds = c(.3,.7)) %>%
filter(type=="After Melding") %>%
mutate(name = gsub("$P(S_1|untested)$",
"$Pr(S_1|untested)$", name,fixed=TRUE),
name = gsub("$P(S_0|test+,untested)$",
"$Pr(S_0|test+,untested)$", name, fixed=TRUE)) %>%
mutate(name = factor(name,
levels = c(
"$\\alpha$",
"$\\beta$",
"$Pr(S_1|untested)$",
"$Pr(S_0|test+,untested)$")))
plt1 <- theta %>%
filter(name != "$M(\\theta) = Pr(S_0|test+,untested)$" & name != "$Pr(S_0|test+,untested)$") %>%
ggplot() +
geom_ribbon(aes(x=x, ymin=0, ymax= value, fill='Before Melding'),alpha=.7, color=NA)+
geom_density(aes(x=value, fill = 'After Melding'),
color=NA,
data = melded_long[melded_long$name !=
"$Pr(S_0|test+,untested)$",]) +
geom_point(aes(x=1,y=0),size=0, color=NA) +
facet_wrap(~name,
labeller= as_labeller(TeX, default = label_parsed),
ncol = 3, scales="free") +
theme_c()+
theme(legend.position="right",
axis.text.x=element_text(size=6),
axis.title.x=element_text(size=12),
axis.title.y=element_text(size=12)) +
theme(legend.position="none") +
scale_x_continuous(n.breaks=5) +
labs(x="Value")+
scale_fill_manual(values = c("Before Melding"= "#5670BF", "Induced" = "#B28542",
"After Melding" = "#418F6A"))
theta_samp <- theta_samp %>% mutate(name ="$Pr(S_0|test+,untested)$" )
plt2 <- theta %>%
filter(name== "$Pr(S_0|test+,untested)$") %>%
ggplot() +
geom_ribbon(aes(x=x, ymin=0, ymax= value, fill='Before Melding'),alpha=.7)+
geom_density(aes(x=value, fill = 'After Melding'),
alpha=.7,color=NA,
data = melded_long[melded_long$name ==
"$Pr(S_0|test+,untested)$",]) +
geom_density(aes(x= phi_induced, fill="Induced"),
color = NA,
data = theta_samp,
alpha =.7) +
facet_wrap(~name,
labeller= as_labeller(TeX, default = label_parsed),
ncol = 3, scales="free") +
theme_c() +
theme(legend.position="right",
axis.text.x=element_text(size=6),
axis.title.x=element_text(size=12),
axis.title.y=element_text(size=12))+
scale_x_continuous(n.breaks=5, limits =c(0,1)) +
labs(x="Value",
fill = "")+
scale_fill_manual(values = c("Before Melding"= "#5670BF",
"Induced" = "#B28542",
"After Melding" = "#418F6A"))
plt1 / plt2